In [1]:
import string
import scipy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
import sys
import glob
import cPickle
import skimage
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
import common_functions as cf
import pims
import passing_event_functions as pef
import trackpy_helper_functions as tphf
import optical_binding_of_driven_particles as obdp
import half_nanoplate_functions as hnf
Load the database
In [2]:
store = pd.HDFStore("K:\Pat's_Projects\ParticleTrajectoryData\passing_event_data_trackpy_refined_extra.h5", mode='r')
keys = store.index['key']
Dwell Times with Filtering Particles Aligned on Ring Trap¶
In [5]:
drop_cols = ['x pos', 'y pos', 'nn_num', 'nn_dist', 'nn_id', 'theta_nn_num', 'theta_nn_dist']
In [6]:
all_l_dwells = []
Ls = [0, 1, 2, 3, 4, 5, -5]
for L in Ls:
keys = store.index[store.index.L == L]['key']
one_l_dwell_time = []
print L
for idx, key in keys.iteritems():
df = store.get(key)
valid_region = df.query('-15 < del_theta < 15')
# valid_region = valid_region[(abs(valid_region.dist_avg_r) < 3) | (df.last_pass_event == True)]
# valid_region = valid_region[(abs(valid_region.dist_avg_r_other_part) < 3) | (df.last_pass_event == True)]
pass_success = df[df.last_pass_event==True].drop_duplicates(['frame', 'pair_id'])
pass_success = pass_success.query('-40 < del_r < 40')
target_indicies = pass_success.index
#print target_indicies.sort_values(),key
for target in target_indicies:
success_pass = df.loc[target]
pass_frame = success_pass.frame
track_id = success_pass['track id']
pair_id = success_pass.pair_id
valid_frames = df[df.frame <= pass_frame]
valid_frames = valid_frames[valid_frames['track id'] == track_id]
valid_frames = valid_frames[valid_frames.pair_id == pair_id]
valid_frames = valid_frames[valid_frames.pass_event_start == True]
start_frame = valid_frames.sort_values('frame').iloc[-1]
start_index = valid_frames.sort_values('frame').index[-1]
frames_dwell = df[df.frame <= start_frame.frame]
frames_dwell = frames_dwell[((frames_dwell.del_theta > -15) & (frames_dwell.del_theta < 15)) | (frames_dwell.pass_event_start==True)]
start_to_pass = pass_frame - start_frame.frame
#print frames_dwell, start_index
dwell_time = hnf.trajectory_dwell_time_to_index(frames_dwell, start_index, max_frames_absent=3)
one_l_dwell_time.append(dwell_time+start_to_pass)
all_l_dwells.append(one_l_dwell_time)
In [7]:
single_exp = lambda x, p1,p2: p1*np.exp(-p2*x)
Ls = [0, 1, 2, 3, 4, 5, -5]
for num, data in enumerate(all_l_dwells):
if len(data) == 0:
continue
print 'L='+str(Ls[num])
bin_num = hnf.hist_bin_optimization(data)
plt.figure()
plt.hist(data, bins=70)
plt.show()
# Histogram
bins = hnf.hist_bin_optimization(data)
hist, bins = np.histogram(data, bins=70)
mid_bins = (bins[:-1] + bins[1:])/2.0
d_t = 1 - np.cumsum(hist)/float(sum(hist))
params, cov = scipy.optimize.curve_fit(single_exp, mid_bins, d_t, p0=[0.5,0.0095], maxfev=2000)
# else:
# params, cov = scipy.optimize.curve_fit(single_exp, mid_bins, d_t, maxfev=2000)
fit_x = np.linspace(min(mid_bins)*0.90, max(mid_bins)*1.10, 300)
fit_y = [single_exp(xp, params[0], params[1]) for xp in fit_x]
# Plot
fig = plt.figure(figsize=(3.5,2.5))
ax = plt.subplot()
plt.plot(mid_bins/110, d_t, 'o')
plt.plot(fit_x/110, fit_y, 'k')
#plt.text(0.7, 0.7, str(params[1]*(1/110.0)), transform=ax.transAxes)
#plt.title(Ls[num])
#plt.title('Cumulative Dwell Time L='+str(Ls[num]))
plt.xlabel('Dwell Time (sec)')
plt.ylabel('$1-\mathregular{CDF}$')
#fig.savefig('passing_dwell_time_L'+str(num)+'.png', bbox_inches='tight', dpi=800)
plt.show()
In [8]:
f = open("C:\Users\Scherer Lab E\Box Sync\Passing\Figures and Data\Data\pass_event_dwell_time.pkl", mode='w')
cPickle.dump(all_l_dwells, f)
f.close()
Number of Passing Events Per L¶
In [17]:
index_df = store.index.copy()
for idx, series in index_df.iterrows():
key = series['key']
df = store.get(key)
num_pass_events = len(df[df.last_pass_event == True])
index_df.loc[idx, 'num_events'] = num_pass_events
index_df.loc[idx, 'num_events_per_frame'] = num_pass_events/float(df.frame.max())
num_particles = int(np.round(series.avg_particles))
combinations = scipy.misc.comb(num_particles, 2)
index_df.loc[idx, 'comb_particles'] = combinations
index_df.loc[idx, 'true_num_frames'] = len(df.groupby('frame'))
index_df.loc[idx, 'true_events_per_frame'] = num_pass_events/float(len(df.groupby('frame')))
normed_num_events = index_df.num_events / index_df.comb_particles
index_df['normed_num_events'] = normed_num_events
normed_num_events = index_df.num_events_per_frame / index_df.comb_particles
index_df['normed_num_events_per_frame'] = normed_num_events
true_normed_events = index_df.true_events_per_frame / index_df.avg_particles
index_df['true_normed_events_per_frame'] = true_normed_events
Number of Passing Events with L¶
In [14]:
plt.plot(index_df['L'], index_df['num_events'], 'o')
plt.ylabel('Number of Passing Events')
plt.xlabel('Angular number, L')
Out[14]:
Normalized Number of Passing Events with L¶
The number of passing events is divided by the combination of the number of possible pairs of particles in the experiment.
In [15]:
plt.plot(index_df['L'], index_df['normed_num_events'], 'o')
plt.ylabel('Normalized Number of Passing Events')
plt.xlabel('Angular number, L')
Out[15]:
Average Normalized Number of Passing Events with L¶
In [16]:
grouped = index_df.groupby('L')
final_l = []
final_num_pass = []
for name, group in grouped:
final_l.append(name)
finite_num_events = group[group.normed_num_events < np.inf]
final_num_pass.append(finite_num_events.normed_num_events.mean())
plt.plot(final_l, final_num_pass, 'o')
plt.ylabel('Avg Normalized Number of Passing Events')
plt.xlabel('Angular number, L')
Out[16]:
Number of Passing Events per Frame with L¶
In [17]:
plt.plot(index_df['L'], index_df['num_events_per_frame'], 'o')
plt.ylabel('Number of Passing Events')
plt.xlabel('Angular number, L')
Out[17]:
Normalized Number of Passing Events per Frame with L¶
The number of passing events is divided by the combination of the number of possible pairs of particles in the experiment.
In [18]:
plt.plot(index_df['L'], index_df['normed_num_events_per_frame'], 'o')
plt.ylabel('Normalized Number of Passing Events per Frame')
plt.xlabel('Angular number, L')
Out[18]:
In [15]:
f = open("C:\Users\Scherer Lab E\Box Sync\Passing\Figures and Data\Data\pass_events_per_sec_normed.pkl", mode='w')
cPickle.dump([index_df['L'], index_df['normed_num_events_per_frame']], f)
f.close()
Average Normalized Number of Passing Events with L¶
In [16]:
grouped = index_df.groupby('L')
final_l = []
final_num_pass = []
for name, group in grouped:
final_l.append(name)
finite_num_events = group[group.normed_num_events < np.inf]
final_num_pass.append(finite_num_events.normed_num_events_per_frame.mean())
plt.plot(final_l, final_num_pass, 'o')
plt.ylabel('Avg Normalized Number of Passing Events')
plt.xlabel('Angular number, L')
Out[16]:
Normalized (by number of particles) Number of Passing Events per True Frames with L¶
The number of passing events is divided by the number of particles in the experiment.
In [22]:
plt.plot(index_df['L'], index_df['true_events_per_frame']/ index_df.comb_particles,'o')
plt.ylabel('Normalized Number of Passing Events per True Frames')
plt.xlabel('Angular number, L')
Out[22]:
Normalized (by number of particles) Number of Passing Events per True Frames with L¶
The number of passing events is divided by the number of particles in the experiment.
In [23]:
plt.plot(index_df['L'], index_df['true_normed_events_per_frame'], 'o')
plt.ylabel('Normalized Number of Passing Events per Frame')
plt.xlabel('Angular number, L')
Out[23]:
Average Normalized (by number of particles) Number of Passing Events with L¶
In [20]:
grouped = index_df.groupby('L')
final_l = []
final_num_pass = []
for name, group in grouped:
final_l.append(name)
finite_num_events = group[group.true_normed_events_per_frame < np.inf]
final_num_pass.append(finite_num_events.true_normed_events_per_frame.mean())
plt.plot(final_l, final_num_pass, 'o')
plt.ylabel('Avg Normalized Number of Passing Events')
plt.xlabel('Angular number, L')
Out[20]:
Look at travel distances making sure the particles are near the center of the trap when considered entering¶
In [27]:
def distance_traveled_before_passing_within_separation_near_ring_center(df, del_theta_separation, min_r=50, deviation_from_cent=2):
df_copy = df.copy()
df_copy['begin_paired_up_frame'] = np.nan
df_copy['distance_traveled_as_pair'] = np.nan
df_events = df[(df.last_pass_event == True) & (df.del_theta < 0) & (abs(df.del_r) < min_r)]
for idx, series in df_events.iterrows():
track_id = series['track id']
pair_id = series['pair_id']
last_pass_frame = series['frame']
df_valid = df_copy.query('frame < @last_pass_frame')
df_valid = df_valid[df_valid.pair_id == pair_id]
df_valid = df_valid[df_valid['track id'] == track_id]
#print df_valid.query('del_theta < @del_theta_separation').drop(['frame', 'track id', 'x pos', 'y pos', 'nn_num', 'nn_id', 'nn_dist'], axis=1)
df_less_than_sep = df_valid.del_theta < -del_theta_separation
df_greater_than_sep = df_valid.del_theta > -del_theta_separation
df_particle_close_trap = abs(df_valid.dist_avg_r) < deviation_from_cent
df_other_particle_close_trap = abs(df_valid.dist_avg_r_other_part) < deviation_from_cent
#print df_less_than_sep, df_greater_than_sep
entering_range = df_less_than_sep.shift(1) * df_greater_than_sep * df_particle_close_trap * df_other_particle_close_trap
if len(entering_range[entering_range == 1]) == 0:
continue
latest_entering_range_index = entering_range[entering_range == 1].index[-1]
#print latest_entering_range_index
frame_begin_pair = df_valid.loc[latest_entering_range_index, 'frame']
assignment_sel = (df.frame == last_pass_frame) & (df.pair_id == pair_id)
df_copy.loc[assignment_sel, 'begin_paired_up_frame'] = frame_begin_pair
df_find_travel = df.query('@frame_begin_pair <= frame <= @last_pass_frame')
df_find_travel = df_find_travel[df_find_travel.pair_id == pair_id]
df_find_travel = df_find_travel[df_find_travel['track id'] == track_id]
theta_midpoint = ((df_find_travel.del_theta/2.0) + df_find_travel.theta) % 360
theta_disp = theta_midpoint.shift(-1) - theta_midpoint
theta_disp = theta_disp - 360 * np.round(theta_disp/360)
df_copy.loc[assignment_sel, 'distance_traveled_as_pair'] = theta_disp.sum()
return df_copy
In [28]:
all_l_pass_dist = []
Ls = [0, 1, 2, 3, 4, 5, -5]
for L in Ls:
keys = store.index[store.index.L == L]['key']
one_l_pass_dist = pd.Series()
print L
for idx, key in keys.iteritems():
df = store.get(key)
df_new = distance_traveled_before_passing_within_separation_near_ring_center(df, 10)
df_new_dropped = df_new.drop_duplicates(['distance_traveled_as_pair'])
travel_distance = df_new_dropped['distance_traveled_as_pair'].dropna()
one_l_pass_dist = pd.concat([one_l_pass_dist, travel_distance])
all_l_pass_dist.append(one_l_pass_dist)
In [29]:
single_exp = lambda x, p1,p2: p1*np.exp(-p2*x)
Ls = [0, 1, 2, 3, 4, 5, -5]
for num, data in enumerate(all_l_pass_dist):
if len(data) == 0:
continue
print 'L='+str(Ls[num])
bin_num = hnf.hist_bin_optimization(data)
plt.figure()
if num == 6:
plt.hist(data, bins=70, range=[-200,0])
print data.mean()
else:
plt.hist(data, bins=70, range=[0,200])
print data.mean()
plt.xlabel('Pair Travel Distance (degrees)')
plt.ylabel('Counts')
plt.show()
# # Histogram
# bins = hnf.hist_bin_optimization(data)
# hist, bins = np.histogram(data, bins=70)
# mid_bins = (bins[:-1] + bins[1:])/2.0
# d_t = 1 - np.cumsum(hist)/float(sum(hist))
# params, cov = scipy.optimize.curve_fit(single_exp, mid_bins, d_t, p0=[0.5,0.0095], maxfev=2000)
# # else:
# # params, cov = scipy.optimize.curve_fit(single_exp, mid_bins, d_t, maxfev=2000)
# fit_x = np.linspace(min(mid_bins)*0.90, max(mid_bins)*1.10, 300)
# fit_y = [single_exp(xp, params[0], params[1]) for xp in fit_x]
# # Plot
# plt.figure()
# ax = plt.subplot()
# plt.plot(mid_bins, d_t, 'o')
# plt.plot(fit_x, fit_y, 'k')
# plt.text(0.7, 0.7, str(params[1]), transform=ax.transAxes)
# #plt.title(Ls[num])
# #plt.title('Cumulative Dwell Time L='+str(Ls[num]))
# plt.xlabel('Dwell Time (sec)')
# plt.ylabel('$1-\mathregular{CDF}$')
# plt.show()
